An immune-related microRNA signature prognostic model for pancreatic carcinoma and association with immune microenvironment

To establish a prognostic model based on immune-related microRNA (miRNA) for pancreatic carcinoma. Weighted correlation network analysis (WGCNA) was performed using the "WGCNA" package to find the key module genes involved in pancreatic carcinoma. Spearman correlation analysis was conducted to screen immune-related miRNAs. Uni- and multi-variate COX regression analyses were carried out to identify miRNAs prognostic for overall survival (OS) of pancreatic carcinoma, which were then combined to generate a prognostic model. Kaplan–Meier survival analysis, receiver operating characteristic (ROC) analysis, distribution plot of survival status in patients and regression analysis were collectively performed to study the accuracy of the model in prognosis. Target genes of the miRNAs in the model were intersected with the key module genes, and a miRNA–mRNA network was generated and visualized by Cytoscape3.8.0. TIMER analysis was conducted to study the abundance of immune infiltrates in tumor microenvironment of pancreatic carcinoma. Expression levels of immune checkpoint genes in subgroups stratified by the model were compared by Wilcoxon test. Gene Set Enrichment Analysis (GSEA) was performed to analyze the enriched signaling pathways between subgroups. Differential analysis revealed 1826 genes differentially up-regulated in pancreatic carcinoma and 1276 genes differentially down-regulated. A total of 700 immune-related miRNAs were obtained, of which 7 miRNAs were significantly associated with OS of patients and used to establish a prognostic model with accurate predictive performance. There were 99 mRNAs overlapped from the 318 target genes of the 7 miRNAs and the key modules genes analyzed by WGCNA. Patient samples were categorized as high or low risk according to the prognostic model, which were significantly associated with dendritic cell infiltration and expression of immune checkpoint genes (TNFSF9, TNFRSF9, KIR3DL1, HAVCR2, CD276 and CD80). GSEA showed remarkably enriched signaling pathways in the two subgroups. This study identified an immune-related 7-miRNA based prognostic model for pancreatic carcinoma, which could be used as a reliable tool for prognosis.

Pancreatic carcinoma (PC) is one of the most fatal malignancies associated with a 5-year survival rate of 8% 1 . First presenting symptoms typically include abdominal or back pain, obstructive jaundice and weight loss. Currently, PC is mainly diagnosed by computed tomography (CT) and treated via surgery. It is poorly sensitive to chemotherapeutic agents. From early asymptomatic in most cases until tumor invasion to surrounding tissues or distance metastasis, approximately 80-85% of the patients have an advanced, un-resectable tumor at the time of diagnosis 2 . In this context, early diagnosis and accurate prognosis are of vital significance for effective management of PC.
MicroRNA (miRNA) is a class of small, non-coding RNA molecules with a size of 17-25 nucleotides 3 . It was first reported by Ambros and Ruvkun early in 1993 4 . It is well believed that miRNA-induced silencing www.nature.com/scientificreports/ age linkage hierarchical clustering was performed based on the differential TOM measures to classify genes of similar expression patterns into a gene module. P value in t test was calculated and visualized.
Screening of immune-related miRNAs prognosis for PC. Single-sample gene set enrichment analysis (ssGSEA) was performed in samples of the TCGA-PAAD dataset to obtain immune scores of M13664 (immune system process) and M19817 (immune response) gene sets based on the MSigDB database (http:// softw are. broad insti tute. org/ gsea/ index. jsp) [15][16][17] . Spearman correlation coefficient was calculated to screen immune scorerelated miRNA (|R|> 0.3, P < 0.05), which were then intersected with the prognostic miRNAs analyzed by univariate analysis to obtain prognostic immune-related miRNAs.
Construction and validation of immune-related miRNA prognostic model. Immune-related miRNAs prognosis for OS of PC were screened by uni-variate analysis. Further multi-variate COX regression analysis was performed to identify signature miRNAs and construct a miRNA prognostic model(R package "survival") 18 . A risk scoring system was established based on the model and defined as a sum of the product of expression level of each miRNA and corresponding regression coefficient. Each patient was scored and divided into two cohorts by taking the median risk score as the threshold. KM survival analysis was applied to study the prognostic significance of the risk score 19 . ROC curve and distribution plot of survival status in the two cohorts were generated to assess the accuracy of the risk score in prognosis(R package "survivalROC") 20 . Independence of the risk score in prognosis was evaluated using uni-and multi-variate analyses, with the involvement of some other significant clinical features.

Establishment of immune-related miRNA-mRNA interaction network.
Target genes for the model miRNAs were predicted using the miRDB, miRTarBase and Target Scan databases respectively [21][22][23] . The predicted genes were then intersected with the genes in the key module identified by WGCNA. A miRNA-mRNA interaction network was then generated and visualized by Cytoscape3.8.0. Different regulatory relationships were distinguished by colors 24 .
Immune cell infiltration and immune checkpoint. CIBERSORT is a deconvolution algorithm developed by Bindea G et al. It can estimate the cell composition and quantify the abundance of specific cell groups in a complex environment using standardized gene expression data 25 . Here, transcriptome data derived from the TCGA-PAAD database were sorted and normalized. Relative abundance of infiltration of 22 immune cells in TME was estimated based on the normalized data using TIMER to study the relationship between the risk score and immune cell infiltration. Since tumor immune therapy has become more popular, expression of immune checkpoint genes in two groups was compared to study the sensitivity of patients to immune therapy.
Statistical analysis. R 4.1.2 and Strawberry Perl fulfilled all statistical analyses of the study. R package "WGCNA" was used to find key module of PC. Spearman correlation analysis was conducted to screen immunerelated miRNAs. The model prognostic was established based on multivariate Cox proportional hazards regression model. Both the univariate and the multivariate Cox regression analysis were performed with the "survival" package. The ROC curves were drawn by the package of "survivalROC" in R. Cytoscape3.8.0 software was applied to generate miRNA-mRNA interaction network. TIMER analysis was carried out to study the abundance of immune infiltrates in TME of PC. Expression levels of immune checkpoint genes in subgroups stratified by the model were compared by Wilcoxon test. GSEA was performed to analyze the enriched signaling pathways between subgroups. P < 0.05 was considered statistically significant.

Results
Key module identified by WGCNA. Following WGCNA, 13 gene modules were initially obtained. Further splicing based on Pearson correlation analysis demonstrated 10 gene modules eventually, which were displayed in a Heat map (Fig. 2). The Brown module (r = 0.2, P = 0.07) and Green module (r = − 0.26, P = 4e − 04) were found to be highly associated with PC. Specifically, the Brown module containing 1,826 genes was positively associated with PC, while the Green module composed of 1276 genes was negatively associated with PC. Genes of the two modules were selected for further analysis.
Screening of immune-related miRNAs. Spearman correlation coefficient was calculated, and a total of 700 immune-related miRNAs were obtained.
Construction and validation of immune-related miRNA prognostic model. 700 immune-related miRNAs were randomly divided into train and test sets. Uni-variate COX regression analysis was performed to obtain 20 candidate immune-related miRNAs of prognostic significance in PC (Table 1). Subsequently, multi-variate COX regression models were established and a 7-miRNA based prognostic model was identified (Table 2). A risk score was accordingly established and formulated as miR-550a-3-5p*-0.472400389 + miR-  www.nature.com/scientificreports/ was found that all these 7 miRNAs were closely associated with the OS of patients, including miR-550a-3-5p, miR-3613-5p, miR-491-3p and miR-1179 showing a negative correlation and miR-221-3p, miR-424-5p and miR-3614-3p showing a positive correlation (Fig. 3). KM survival analysis was also performed in the high-and low-risk cohorts of the train and test sets respectively. Significant difference in OS of the two groups in both sets was observed (P < 0.001) (Fig. 4). AUC values of the ROC curves for survival were 0.715, 0.754 and 0.678 respectively, demonstrating high accuracy of the performance of the model (Fig. 5). Additionally, distribution plot of the survival status of patients in the two groups revealed an increased number of death with the increase of risk score (Fig. 6). Finally, the risk score was detected to be prognosis for OS of PC as detected by uni-variate analysis (P < 0.001). According to P < 0.2, further multi-variate analysis indicated that the risk score was independent of other factors in prognosis of PC (P < 0.001) (Fig. 7).

Establishment of immune-related miRNA-mRNA interaction network.
Target genes for the model miRNAs were predicted using the miRDB, miRTarBase and TargetScan databases respectively. The intersected genes were selected for further analysis, including 100 mRNAs for miR-221-3p, 272 mRNAs for miR-424-5p, 15 mRNAs for miR-491-3p, 14 mRNAs for miR-550a-3-5p, 11 mRNAs for miR-1179, 11 mRNAs for miR-3613-5p, and 7 mRNAs for miR-3614-3p (Fig. 8). After intersection with the key module genes in WGCNA, a total of 99 genes were obtained and considered significantly associated with the immunity and prognosis of PC (Fig. 9). A miRNA-mRNA interaction network was constructed based on the 7 model miRNAs and 99 mRNAs, which the red represented a positive correlation and the blue represented a negative correlation (Fig. 10).
Immune cell infiltration and immune checkpoint. Composition of immune cell subgroups in TME of the high-and low-risk groups was analyzed and the infiltration abundance was displayed by a Heat map  www.nature.com/scientificreports/ (Fig. 11). Infiltration of dendritic cells was demonstrated important between the two groups. Moreover, expression of immune checkpoint genes was examined. It was found that the risk score was associated with expression levels of TNFSF9, TNFRSF9, KIR3DL1, HAVCR2, CD276, especially CD80 (Fig. 12).
GSEA for target genes of the model miRNAs. By GSEA, the gene set at high risk was mainly involved in the pathways of cell cycle, chronic myeloid leukemia, prostate cancer, small cell lung cancer,steroid hormone biosynthesis. In the meantime, the gene set at low risk was predominantly associated with pathways involved in   www.nature.com/scientificreports/ steroid hormone biosynthesis,calcium signaling pathway,long term potentiation,maturity onset diabetes of the young,neuroactive ligand receptor interaction (Fig. 13).

Discussion
This study, for the first time, constructed a prognostic model using immune-related miRNAs and validated its accurate, potent performance. We firstly analyzed the key modules associated with prognosis of PC using WGCNA, and then adopted differential analysis to obtain 1826 up-regulated and 1276 down-regulated genes in PC comparing with normal tissues. Spearman correlation coefficient was calculated and 700 immune-related miRNAs were found. Uni-and multi-variate COX regression analyses were performed to screen immune-related miRNA associated with OS. Eventually, a 7-miRNA-based prognostic model for PC was generated and a modelbased risk scoring system was established with excellent performance detected by KM, ROC, distribution plot of survival status in subgroups, and uni-and multi-variate COX regression analyses. Furthermore, a miRNA-mRNA network was built using the Cytoscape3.8.0 software, according to the 99 target genes obtained from the predicted mRNAs of the 7 model miRNAs and key module genes. Moreover, immune cell infiltration was analyzed using the TIMER in the subgroups stratified by the model-based risk score, and dendritic cells were found as the most significantly different immune infiltrate in TME of the subgroups. Wilcoxon test also indicated correlation of the risk score with expression of immune checkpoint genes, including TNFSF9, TNFRSF9, KIR3DL1, HAVCR2, CD276, especially CD80. GSEA showed the remarkably enriched signaling pathways in the two subgroups.
There have been multiple studies devoted to constructing prognostic models for PC. For instance, Weng W et al. 26 established a prognostic model for PC using 14 significant mRNAs; Yan et al. 27 reported a 4-gene-based  Currently, prognostic miRNAs for PC have been extensively studied. Guo et al. 31 reported that high expression levels of miR-21, miR-212, miR-675, miR-142-5p, miR-204 and low expression levels of miR-148a, miR-187, let-7g were associated with poor prognosis of PC. In that study, they also found several other miRNAs of prognostic significance, including miR-30a-3p, miR-105, miR-127, miR-187, miR-452, miR-518a-2, miR-155,   www.nature.com/scientificreports/ increased expression of miR-21 and miR155 was associated with the decline in survival, liver metastasis, lymph node status and increased resistance to Gemcitabine.
In the current study, we successfully established a prognostic model using 7 immune-related miRNAs. Lu et al. 34 found that expression of miR-550a in severe acute pancreatitis was down-regulated in patients combining with acute lung injury compared with patients without acute lung injury. Qin et al. 35 revealed that miR-3613-5p was prognostic for renal clear cell carcinoma. Ma et al. 36 also reported miR-3613-5p as a prognostic biomarker for pancreatic carcinoma and found that the target genes of miR-3613-5p might be correlated with the p53 signaling pathway. Research also found that miR-3613-5p was a key miRNA in cell malignancy of liver cancer induced by aflatoxin B1 37 . Wang et al. 38 found that miR-221-3p was independent of the prognosis of hepatocellular carcinoma (HCC), while Wang et al. 39 reported the key role of miR-221-3p in thyroid cancer. Moreover, Fang et al. 40 established a prognostic model for breast cancer based on 13 miRNAs including miR-221-3p. Abak et al. 41 analyzed breast biopsy samples and found increased expression of miR-221-3p in tumor tissue compared to that in adjacent normal tissue. Xie et al. 42 adopted high-throughput sequencing from peripheral blood mononuclear cells in small cell osteosarcoma and detected dysregulation of miR-221-3p. Kandhavelu et al. 43 found that miR-424-5p was associated with 10 biomarkers of colon carcinoma, while only two of them, including microtubule-associated protein-2 (MAP2) and cyclin gene (CCN) D1, were experimentally validated. Liang et al. 44 identified miR-424-5p as a potential prognostic biomarker for gastric cancer. In the meantime, Liu et al. 45 experimentally validated that the up-regulation of miR-424-5p induced by astragaloside IV could inhibit the epithelial-mesenchymal transition (EMT) and angiogenesis in gastric cancer, thereby to play a role in treatment of this cancer. Wang et al. 46 found and validated the role of miR-424-5p as a factor in prognosis of HCC. Xu et al. 47 also reported miR-424-5p played a vital role in onset of bile duct carcinoma. Additionally, miR-424-5p might reverse the progression of thyroid cancer via regulating clusterin (CLU) and apolipoprotein (APO) 48 . Ranjha et al. 49 found decreased expression of miR-491-3p in ulcerative colitis located in the rectal sigmoid colon region than expression in ulcerative colitis located in the ascending colon. circANKRD36 might interact with miR-3614-3p to participate in Type 2 Diabetes and inflammation, supporting their role as potential biomarkers 50 . miR-1179 inhibited the in vivo growth of PC by down-regulating the expression of E2F transcription factor 5 51 . Moreover, miR-1179 could suppress the growth and invasion of thyroid, gastric, esophageal squamous-cell carcinoma, non-small cell lung, cervical, breast, and nasopharyngeal carcinomas via regulating downstream target genes [52][53][54][55][56][57][58] .
To conclude, this is the first study that constructed a prognostic model for PC using immune-related miRNAs, aiming to provide a new direction for clinical prognosis of PC. The current study still has some deficiencies, such as lack of validation by molecular biology or relevant clinical trials. Additionally, the miR-491-3p and miR-550a-3-5p have not been fully studied in the field of tumor, requiring further in-depth research.

Data availability
TCGA-PAAD datasets is downloaded from the The Cancer Genome Atlas Program (TCGA) database (https:// portal. gdc. cancer. gov/ proje cts/ TCGA-PAAD). All data and R script in this study are available from the corresponding author on reasonable request.